####################################################
####          Replication Code for              ####
####     Bayesian versus maximum likelihood     ####
####      estimation of treatment effects       ####
####          in bivariate probit               ####  
####       instrumentalvariable models          ####    
####################################################
####                                            ####  
####                3/ 7/2017                   ####         
####     this file runs rcode to create Figures ####  
####            in Appendix F                   ####  
####################################################




#### recreate simulation parameters
N_seq = c(50, 100, 250, 500)
cor = c(0.25,0.5,0.75)
IV_strength = c(0.5, 1, 1.5, 2)
treatment = c(0, 1, 2, 3, 4)
simulation_para = data.frame(do.call(rbind,replicate(500,as.matrix(expand.grid(N_seq, cor, IV_strength, treatment)),simplify = F)))
names(simulation_para) = c("N","cor","IVcoef","treatment")
simulation_para = simulation_para[order(simulation_para$N, simulation_para$cor, simulation_para$IVcoef, simulation_para$treatment),]
simulation_para$iteration = seq(1,120000)



### load results from Bayesian models with different priors

### load Bayesian model results  and Stata output
load("R_output15.rda")
load("R_output25.rda")
load("R_output35.rda")
load("R_outputUniform.rda")

load("stata_output.rda")
### stata bootstrapped CIs for ATE
load("ate_estimates.rda")
names(ate_est)[1] = "iteration"




##combine stata and bayes results by iteration
stata = data.frame(stata_output[,c("iteration", "ATE")])
names(stata) = c("iteration", "StataAte")
stata = join(stata, ate_est[, c("iteration", "Lower","Upper")], by = "iteration")
names(stata) <- c("iteration", "StataAte", "StataLow95", "StataHigh95")

### bayes res
Bayes_res.25 = data.frame(t(R_output25["ate",c("mean","median", "2.5%","97.5%", "TRUE"),]), (R_output25["iteration","Rhat/Misc",]))
names(Bayes_res.25) = c("BayesMean.25","BayesMedian.25", "BayesLow95.25","BayesHigh95.25", "TrueAte","iteration")

Bayes_res.15 = data.frame(t(R_output15["ate",c("mean","median", "2.5%","97.5%"),]), (R_output25["iteration","Rhat/Misc",]))
names(Bayes_res.15) = c("BayesMean.15","BayesMedian.15", "BayesLow95.15","BayesHigh95.15","iteration")

Bayes_res.35 = data.frame(t(R_output35["ate",c("mean","median", "2.5%","97.5%"),]), (R_output25["iteration","Rhat/Misc",]))
names(Bayes_res.35) = c("BayesMean.35","BayesMedian.35", "BayesLow95.35","BayesHigh95.35","iteration")

Bayes_res.Uniform = data.frame(t(R_outputUniform["ate",c("mean","median", "2.5%","97.5%"),]), (R_output25["iteration","Rhat/Misc",]))
names(Bayes_res.Uniform) = c("BayesMean.Uniform","BayesMedian.Uniform", "BayesLow95.Uniform","BayesHigh95.Uniform","iteration")






sim_results = join(Bayes_res.25, Bayes_res.15, by = "iteration")
sim_results = join(sim_results, Bayes_res.35, by = "iteration")
sim_results = join(sim_results, Bayes_res.Uniform, by = "iteration")
sim_results = join(sim_results, stata, by = "iteration")
sim_results = join(sim_results, simulation_para, by = "iteration")





### create indicator whether 95CI includes true value
### create indicator whether 95CI includes true value
sim_results$Bayes95Contain.15 = ifelse(sim_results$TrueAte <= sim_results$BayesHigh95.15 & sim_results$TrueAte >= sim_results$BayesLow95.15, 1, 0)
sim_results$Bayes95Contain.25 = ifelse(sim_results$TrueAte <= sim_results$BayesHigh95.25 & sim_results$TrueAte >= sim_results$BayesLow95.25, 1, 0)
sim_results$Bayes95Contain.35 = ifelse(sim_results$TrueAte <= sim_results$BayesHigh95.35 & sim_results$TrueAte >= sim_results$BayesLow95.35, 1, 0)
sim_results$Bayes95Contain.Uniform = ifelse(sim_results$TrueAte <= sim_results$BayesHigh95.Uniform & sim_results$TrueAte >= sim_results$BayesLow95.Uniform, 1, 0)
sim_results$Stata95Contain = ifelse(sim_results$TrueAte <= sim_results$StataHigh95 & sim_results$TrueAte >= sim_results$StataLow95, 1, 0)


### calculate errors
sim_results$Stata_abs_error = abs(sim_results$StataAte - sim_results$TrueAte)
sim_results$Bayes_abs_error.15 = abs(sim_results$BayesMedian.15 - sim_results$TrueAte)
sim_results$Bayes_abs_error.25 = abs(sim_results$BayesMedian.25 - sim_results$TrueAte)
sim_results$Bayes_abs_error.35 = abs(sim_results$BayesMedian.35 - sim_results$TrueAte)
sim_results$Bayes_abs_error.Uniform = abs(sim_results$BayesMedian.Uniform - sim_results$TrueAte)


#### check for any simulations for stata for which we couldn't calculate the CI
dim(sim_results[is.na(sim_results$StataHigh95) == T & is.na(sim_results$Stata_abs_error)==F,])
sim_results$StataAte = ifelse(is.na(sim_results$StataHigh95) == T, NA, sim_results$StataAte)

### analyze results
#### by N and IV strength
N = c(50, 100, 250, 500)
IVs =  c(0.5, 1, 1.5, 2.0)
Bayes_mse.15 = Bayes_mse.25 = Bayes_mse.35 = Bayes_mse.Uniform = stata_mse = Bayes_mae.15 = Bayes_mae.25 = Bayes_mae.35 = Bayes_mae.Uniform = stata_mae  = Bayes_pct95.15 = Bayes_pct95.25 = Bayes_pct95.35 = Bayes_pct95.Uniform = Stata_pct95 = n_val = iv_val = NULL

for(i in 1: 4){
    for(j in 1:4){
        Bayes_mse.15 = c(Bayes_mse.15, mean_sq(sim_results$BayesMedian.15[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE]))
        Bayes_mse.25 = c(Bayes_mse.25, mean_sq(sim_results$BayesMedian.25[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE]))
        Bayes_mse.35 = c(Bayes_mse.35, mean_sq(sim_results$BayesMedian.35[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE]))
        Bayes_mse.Uniform = c(Bayes_mse.Uniform, mean_sq(sim_results$BayesMedian.Uniform[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE]))
        stata_mse =  c(stata_mse , mean_sq(sim_results$StataAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE]))

        
        Bayes_mae.15 =  c(Bayes_mae.15, MAE(sim_results$BayesMedian.15[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE]))
        Bayes_mae.25 =  c(Bayes_mae.25, MAE(sim_results$BayesMedian.25[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE]))
        Bayes_mae.35 =  c(Bayes_mae.35, MAE(sim_results$BayesMedian.35[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE]))
        Bayes_mae.Uniform =  c(Bayes_mae.Uniform, MAE(sim_results$BayesMedian.Uniform[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE] , sim_results$TrueLate[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE]))
        stata_mae =  c(stata_mae, MAE(sim_results$StataAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE] , sim_results$TrueAte[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE]))

        Bayes_pct95.15 = c(Bayes_pct95.15, mean(sim_results$Bayes95Contain.15[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE]))
        Bayes_pct95.25 = c(Bayes_pct95.25, mean(sim_results$Bayes95Contain.25[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE]))
        Bayes_pct95.35 = c(Bayes_pct95.35, mean(sim_results$Bayes95Contain.35[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE]))
        Bayes_pct95.Uniform = c(Bayes_pct95.Uniform, mean(sim_results$Bayes95Contain.Uniform[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE]))
        Stata_pct95 = c(Stata_pct95, mean(sim_results$Stata95Contain[sim_results$N == N[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE]))
        
        
        n_val = c(n_val, N[i])
        iv_val = c(iv_val, IVs[j])
        }
}       



#### create figures 
data_plot = data.frame(mse = c(Bayes_mse.15, Bayes_mse.25, Bayes_mse.35, Bayes_mse.Uniform, stata_mse), c95 = c(Bayes_pct95.15, Bayes_pct95.25, Bayes_pct95.35, Bayes_pct95.Uniform, Stata_pct95), iv= rep(iv_val,5),n = rep(n_val,5), method = c(rep("Bayes.15",16), rep("Bayes.25",16), rep("Bayes.35",16), rep("Bayes.Uniform",16),rep("Stata",16)))

### Figure F1 in Appendix

p1 = ggplot(data = data_plot[data_plot$n == 50,],aes(y = mse, x = iv))
p1 = p1 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p1 = p1 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p1 = p1 + labs(title = "N = 50", y="MSE", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p1 = p1 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p1 = p1 + theme(legend.position = c(.65, .75), legend.direction = "horizontal", legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 14), legend.title = element_text(size = 14), panel.border = element_rect(colour = "black", fill=NA, size=3.5))# + scale_y_continuous(limits = c(0, 1))
plot(p1)
leg = g_legend(p1)
p1 = p1 + theme(legend.position="none")


p2 = ggplot(data = data_plot[data_plot$n == 100,],aes(y = mse, x = iv))
p2 = p2 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p2 = p2 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p2 = p2 + labs(title = "N = 100", y="MSE", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p2 = p2 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p2 = p2 + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))
plot(p2)


p3 = ggplot(data = data_plot[data_plot$n == 250,],aes(y = mse, x = iv))
p3 = p3 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p3 = p3 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p3 = p3 + labs(title = "N = 250", y="MSE", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p3 = p3 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p3 = p3 + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))
plot(p3)

p4 = ggplot(data = data_plot[data_plot$n == 500,],aes(y = mse, x = iv))
p4 = p4 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p4 = p4 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p4 = p4 + labs(title = "N = 500", y="MSE", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p4 = p4 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p4 = p4 + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))
plot(p4)


pdf("Figure_F1.pdf", width = 12, height = 12 )
grid.newpage()
pushViewport(viewport(layout = grid.layout(4,3, heights = unit(c(11,1, 11, 2), "null"), width = unit(c(11,1, 11), "null"))))   
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))    
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))         
print(p3, vp = viewport(layout.pos.row = 3, layout.pos.col = 1))         
print(p4, vp = viewport(layout.pos.row = 3, layout.pos.col = 3))
leg$vp = viewport(layout.pos.row = 4, layout.pos.col = c(1:3))
grid.draw(leg) 
dev.off()




### now plot 95% coverage
### Figure F2 in Appendix

p1 = ggplot(data = data_plot[data_plot$n == 50,],aes(y = c95, x = iv))
p1 = p1 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p1 = p1 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p1 = p1 + labs(title = "N = 50", y="95% Interval Coverage", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p1 = p1 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p1 = p1 + theme(legend.position = c(.65, .75), legend.direction = "horizontal", legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 14), legend.title = element_text(size = 14), panel.border = element_rect(colour = "black", fill=NA, size=3.5))# + scale_y_continuous(limits = c(0, 1))
p1 = p1 + geom_hline(yintercept = 0.95, colour = I("RED"),  size = 1)
p1 = p1 + geom_hline(yintercept = 1, colour = I("BLUE"), size = 1) + scale_y_continuous(limits = c(0.6, 1.005), breaks = c(0.7,0.8,0.9,0.95,1.0),labels = c("0.7","0.8","0.9","0.95","1.0"))
plot(p1)
leg = g_legend(p1)
p1 = p1 + theme(legend.position="none")


p2 = ggplot(data = data_plot[data_plot$n == 100,],aes(y = c95, x = iv))
p2 = p2 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p2 = p2 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p2 = p2 + labs(title = "N = 100", y="95% Interval Coverage", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p2 = p2 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p2 = p2  + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))# + theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
p2 = p2 + geom_hline(yintercept = 0.95, colour = I("RED"), size = 1)
p2 = p2 + geom_hline(yintercept = 1, colour = I("BLUE"),  size = 1) + scale_y_continuous(limits = c(0.6, 1.005), breaks = c(0.7,0.8,0.9,0.95,1.0),labels = c("0.7","0.8","0.9","0.95","1.0"))
plot(p2)


p3 = ggplot(data = data_plot[data_plot$n == 250,],aes(y = c95, x = iv))
p3 = p3 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p3 = p3 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p3 = p3 + labs(title = "N = 250", y="95% Interval Coverage", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p3 = p3 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p3 = p3  + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))#+ theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
p3 = p3 + geom_hline(yintercept = 0.95, colour = I("RED"), size = 1)
p3 = p3 + geom_hline(yintercept = 1, colour = I("BLUE"), size = 1) + scale_y_continuous(limits = c(0.6, 1.005), breaks = c(0.8,0.9,0.95,1.0),labels = c("0.8","0.9","0.95","1.0"))
plot(p3)

p4 = ggplot(data = data_plot[data_plot$n == 500 ,],aes(y = c95, x = iv))
p4 = p4 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p4 = p4 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p4 = p4 + labs(title = "N = 500", y="95% Interval Coverage", x = "Instrument Strength")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p4 = p4 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p4 = p4 + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))#+ theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
p4 = p4 + geom_hline(yintercept = 0.95, colour = I("RED"), size = 1)
p4 = p4 + geom_hline(yintercept = 1, colour = I("BLUE"), size = 1) + scale_y_continuous(limits = c(0.6, 1.005), breaks = c(0.8,0.9,0.95,1.0),labels = c("0.8","0.9","0.95","1.0"))
plot(p4)

pdf("Figure_F2.pdf", width = 12, height = 12 )
grid.newpage()
pushViewport(viewport(layout = grid.layout(4,3, heights = unit(c(11,1, 11, 2), "null"), width = unit(c(11,1, 11), "null"))))   
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))    
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))         
print(p3, vp = viewport(layout.pos.row = 3, layout.pos.col = 1))         
print(p4, vp = viewport(layout.pos.row = 3, layout.pos.col = 3))
leg$vp = viewport(layout.pos.row = 4, layout.pos.col = c(1:3))
grid.draw(leg) 
dev.off()


##### make plots by treatment for ATE
Treat = c(0, 1, 2, 3, 4)
IVs =  c(0.5, 1, 1.5, 2.0)
Bayes_mse.15 = Bayes_mse.25 = Bayes_mse.35 = Bayes_mse.Uniform = stata_mse = Bayes_mae.15 = Bayes_mae.25 = Bayes_mae.35 = Bayes_mae.Uniform = stata_mae  = Bayes_pct95.15 = Bayes_pct95.25 = Bayes_pct95.35 = Bayes_pct95.Uniform = stata_pct95 = treat_val = iv_val = NULL


for(i in 1: 5){
  for(j in 1:4){
    Bayes_mse.15 = c(Bayes_mse.15, mean_sq(sim_results$BayesMedian.15[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE] , sim_results$TrueAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.15)==FALSE]))
    Bayes_mse.25 = c(Bayes_mse.25, mean_sq(sim_results$BayesMedian.25[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE] , sim_results$TrueAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.25)==FALSE]))
    Bayes_mse.35 = c(Bayes_mse.35, mean_sq(sim_results$BayesMedian.35[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE] , sim_results$TrueAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.35)==FALSE]))
    Bayes_mse.Uniform = c(Bayes_mse.Uniform, mean_sq(sim_results$BayesMedian.Uniform[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE] , sim_results$TrueAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$BayesMedian.Uniform)==FALSE]))
    stata_mse =  c(stata_mse , mean_sq(sim_results$StataAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE] , sim_results$TrueAte[sim_results$treatment == Treat[i] & sim_results$IVcoef == IVs[j] & is.na(sim_results$StataAte)==FALSE]))
    
    treat_val = c(treat_val, Treat[i])
    iv_val = c(iv_val, IVs[j])
  }
}       




data_plot = data.frame(mse = c(Bayes_mse.15, Bayes_mse.25, Bayes_mse.35, Bayes_mse.Uniform, stata_mse), treat = rep(treat_val,5),iv = rep(iv_val,5), method = c(rep("Bayes.15",20), rep("Bayes.25",20), rep("Bayes.35",20), rep("Bayes.Uniform",20), rep("Stata",20)))

p1 = ggplot(data = data_plot[data_plot$iv == 0.5,],aes(y = mse, x = treat))
p1 = p1 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p1 = p1 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p1 = p1 + labs(title = "IV = 0.5", y="MSE", x = "Treatment")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p1 = p1 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p1 = p1 + theme(legend.position = c(.65, .75), legend.direction = "horizontal", legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 14), legend.title = element_text(size = 14), panel.border = element_rect(colour = "black", fill=NA, size=3.5))# + scale_y_continuous(limits = c(0, 1))
plot(p1)
leg = g_legend(p1)
p1 = p1 + theme(legend.position="none")

p2 = ggplot(data = data_plot[data_plot$iv == 1.0,],aes(y = mse, x = treat))
p2 = p2 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p2 = p2 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p2 = p2 + labs(title = "IV = 1", y="MSE", x = "Treatment")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p2 = p2 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p2 = p2  + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))#+ theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
plot(p2)


p3 = ggplot(data = data_plot[data_plot$iv == 1.5,],aes(y = mse, x = treat))
p3 = p3 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p3 = p3 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p3 = p3 + labs(title = "IV = 1.5", y="MSE", x = "Treatment")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p3 = p3 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p3 = p3  + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))#+ theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
plot(p3)

p4 = ggplot(data = data_plot[data_plot$iv == 2.0,],aes(y = mse, x = treat))
p4 = p4 + geom_point(aes(colour = method, shape = method, linetype = method), size = 5) + geom_line(aes(colour = method, linetype = method, shape = method),size = 2.5)
p4 = p4 + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p4 = p4 + labs(title = "IV = 2.0", y="MSE", x = "Treatment")+ theme(axis.text.x=element_text(size=20)) + theme(axis.text.y=element_text(size=20))+ theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20), title = element_text(size = 25))
p4 = p4 + scale_colour_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5", "Bayes, Uniform Prior","Stata MLE"), values = c("#377eb8", "#4daf4a", "#984ea3","#ff7f00", "gray10")) + scale_shape_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(16, 17, 18, 19, 20)) + scale_linetype_manual(name = "Estimation Method", labels = c("Bayes, Prior SD = 1.5", "Bayes, Prior SD = 2.5", "Bayes, Prior SD = 3.5",  "Bayes, Uniform Prior","Stata MLE"), values = c(2, 3, 4, 5, 6))
p4 = p4  + theme(legend.position="none") + theme(panel.border = element_rect(colour = "black", fill=NA, size=3.5))#+ theme(legend.position = c(.8, .2),legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 20), legend.title = element_text(size = 20), panel.border = element_rect(colour = "black", fill=NA, size=2.5)) #+ scale_y_continuous(limits = c(0, 6.5))
plot(p4)

pdf("Figure_F3.pdf", width = 12, height = 12 )
grid.newpage()
pushViewport(viewport(layout = grid.layout(4,3, heights = unit(c(11,1, 11, 2), "null"), width = unit(c(11,1, 11), "null"))))   
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))    
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))         
print(p3, vp = viewport(layout.pos.row = 3, layout.pos.col = 1))         
print(p4, vp = viewport(layout.pos.row = 3, layout.pos.col = 3))
leg$vp = viewport(layout.pos.row = 4, layout.pos.col = c(1:3))
grid.draw(leg) 
dev.off()


